function rmsJitter = jitter2(symbolRate,threshold,hysteresis)
% JITTER2.M - Modified version of jitter1 to acquire data directly from the scope
% calling syntax
%      rmsjitter = jitter2(5000,0,0.1)
% or   rmsjitter = jitter2
% in the latter case default parameter values are used.

% Tektronix, Inc. 2001
% Version 7/18/2001

% use default values if function is called without arguments
if nargin < 1
symbolRate = 5000;
threshold = 0;
hysteresis = 0.1;
end 

% This function calculates the RMS jitter in a waveform.
% Jitter is the difference between the actual time an edge occurs and the
% time where it should have been based on the supplied sample rate.
% The waveform is acquired directly from the scope

strCh = 'ch1';
% change the value below to test for different record lengths
recordLen = 400000;
% Use inside the scope with Tek VISA (board 8, primary address 1)
g = visa('tek','GPIB8::1::INSTR');
% (if running MATLAB on a connected PC, change vendor code and/or GPIB 
% descriptor as necessary e.g. 
% g = visa('ni','GPIB0::1::INSTR');
% g = visa('agilent','GPIB0::1::INSTR');

% Make sure the size of the InputBuffer - in bytes - is sufficient.
set(g,'InputBufferSize',recordLen*2);

fopen(g);
idn = query(g,'*IDN?');
fprintf(g,'HEADER OFF');
fprintf(g,['DATA:SOURCE ' strCh]);
fprintf(g,'DATA:ENCDG SRIBINARY;');
fprintf(g,'DATA:WIDTH 2');
fprintf(g,'ACQUIRE:STATE OFF');
fprintf(g,'ACQUIRE:MODE NORMALSAMPLE');
fprintf(g,'ACQUIRE:STOPAFTER SEQUENCE');
fprintf(g,'ACQUIRE:STATE RUN');

while query(g,'BUSY?','%s','%e'); end;
horizLen = query(g,'HORIZONTAL:RECORD?','%s','%e');

fprintf(g,'DATA:START %d', 1);
fprintf(g,'DATA:STOP %d', recordLen);
fprintf(g,'CURVE?');
dummy_string1 = fscanf(g,'%s',2);
dummy_string2 = fscanf(g,'%s',str2num(dummy_string1(2)));
recordLen2Transfer = min(recordLen,horizLen);
[waveform_raw,count] = fread(g,recordLen2Transfer,'int16');
dummy_string3 = fscanf(g,'%s',1);

% get the sampling interval 
sampleInterval = query(g,'WFMOUTPRE:XINCR?','%s','%e')

% Scale the data 
yunit = query(g,'WFMOUTPRE:YUNIT?');
ymult = query(g,'WFMOUTPRE:YMULT?','%s','%e');
yoff = query(g,'WFMOUTPRE:YOFF?','%s','%e');
yzero = query(g,'WFMOUTPRE:YZERO?','%s','%e');
% scale the data to the correct values
waveform = ymult*(waveform_raw - yoff) - yzero;

% find the edges in the supplied waveform
measuredTime = measureEdgeTiming2(waveform,threshold,hysteresis,sampleInterval);

% preallocate space for the clocks array
clocks=zeros(1,length(measuredTime));

% derive the clocks based on the supplied symbol rate
for index = 2:length(measuredTime);
   clocks(index) = (round(symbolRate * (measuredTime(index) - measuredTime(index - 1)))) + clocks(index-1);
end

% fit the derived clocks and the measured time to a straight line
coef = polyfit(clocks, measuredTime, 1);

% coef(2) is the intercept (a) in the form y = a + bx
% coef(1) is the slope (b) in the form y = a + bx
a = coef(2);
b = coef(1);

measuredAverageSymbolRate = 1/b;
measuredSymbolRateError = (measuredAverageSymbolRate - symbolRate)/symbolRate;

subplot(2,1,1);
plot(waveform);
title(['symbol rate error: ', num2str(measuredSymbolRateError * 100), '%']);
xlabel('samples');
% provide the label in units acquired from the scope
% strtok is needed to remove double quotes
ylabel(['waveform amplitude, ' strtok(yunit,'"')]);

reconstructedTime = a + (clocks .* b);

% jitter is the difference between the measured time and the reconstructed time.
jitter = reconstructedTime - measuredTime;

% see the MATLAB function reference for 'norm'
rmsJitter = norm(jitter)/sqrt(length(jitter));

subplot(2,1,2);
plot(reconstructedTime,jitter);
title(['RMS jitter: ', num2str(rmsJitter*1e6), ' \mus']);
xlabel('time in seconds');
ylabel('jitter in \mus');
% force the x-axis limits to be tight so that both plots line up
set(gca,'XLim',[0 count*sampleInterval])
% close the instrument object
fclose(g);
delete(g);


